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(57) ABSTRACT 

An improved model for simulating atmospheric disturbances 
is disclosed. A scale Kolmogorov spectral may be scaled to 
convert the Kolmogorov spectral into a finite energy von 
Karman spectral and a fractional order pole-zero transfer 
function (TF) may be derived from the von Karman spectral. 
Fractional order atmospheric turbulence may be approxi- 
mated with an integer order pole-zero TF fit, and the approxi- 
mation may be stored in memory. 

20 Claims, 12 Drawing Sheets 


1100 




















FIG. 1 



Frequency, Hz 























Frequency, Hz 



















FIG. 7 













Time, sec 







FIG. 10 



Frequency, Hz 






U.S. Patent 


Jan. 20, 2015 


Sheet 11 of 12 


US 8,938,377 B1 









FIG. 12 


U.S. Patent 


Jan. 20, 2015 


Sheet 12 of 12 


US 8,938,377 B1 












US 8,938,377 B1 


1 

ATMOSPHERIC TURBULENCE MODELING 
FOR AEROSPACE VEHICLES: FRACTIONAL 
ORDER FIT 

CROSS REFERENCE TO RELATED 
APPLICATIONS 

This application is a non-provisional of, and claims priority 
to, U.S. Provisional Patent Application Ser. No. 61/488,388, 
filed May 20, 2011, the Subject matter of which is hereby 
incorporated by reference in its entirety. 

ORIGIN OF THE INVENTION 

The invention described herein was made by an employee 
of the United States Government and may be manufactured 
and used by or for the Government for Government purposes 
without the payment of any royalties thereon or therefore. 

FIELD 

The present invention generally pertains to modeling, and 
more specifically, to modeling atmospheric turbulence for the 
design and safety compliance of aerospace vehicles using a 
fractional order fit. 

BACKGROUND 

Atmospheric turbulence models are often beneficial for the 
design and safety compliance of aerospace vehicles. For 
instance, such models are often used for the design of inlet/ 
engine and flight controls, as well as for studying coupling 
between the propulsion and the vehicle structural dynamics 
for supersonic vehicles. Models based primarily on the Kol- 
mogorov spectrum have been previously utilized to model 
atmospheric turbulence. The Kolmogorov spectrum has an 
energy level that approaches infinity as frequency approaches 
zero. This characteristic makes it difficult to implement con- 
ventional models in the time domain. The Kolmogorov model 
has also been extended in the Tank model to develop a base- 
line of atmospheric turbulence for the High Speed Civil 
Transport (HSCT). The Tank model also covers atmospheric 
acoustic wave disturbance modeling utilizing the von Kar- 
man spectral. 

An approximation of the Kolmogorov model that is com- 
monly used with a finite energy spectrum is the von Karman 
type model. However, simulating the von Karman type mod- 
els in the time domain may be difficult due to the fractional 
order of von Karman type models. 

Some have used Dryden model approximation to the frac- 
tional order von Karman atmospheric model. However, the 
Dryden model is second order compared to the 5/3 fractional 
order of the acoustic velocity atmospheric turbulence spec- 
tral. Thus, the Dryden model underestimates the atmospheric 
disturbance, and does so increasingly with higher frequency. 
As such, the faster an aerospace vehicle travels, the more the 
Dryden model will underestimate the atmospheric distur- 
bance. Accordingly, a more accurate model for simulating 
atmospheric disturbances may be beneficial. 

SUMMARY 

Certain embodiments of the present invention may be 
implemented and provide solutions to the problems and needs 
in the art that have not yet been fully solved by conventional 
atmospheric modeling approaches. For example, certain 
embodiments of the present invention utilize a fractional 


2 

order fit to create accurate models of atmospheric turbulence. 
This modeling methodology may be implemented in soft- 
ware, hardware, or a combination of the two. 

In one embodiment of the present invention, a computer 
5 program embodied on a computer-readable storage medium 
is configured to cause a processor to scale a Kolmogorov 
spectral to convert the Kolmogorov spectral into a finite 
energy von Karman spectral and derive a fractional order 
pole-zero transfer function (TF) from the von Karman spec- 
io tral. The computer program is also configured to cause the 
processor to approximate fractional order atmospheric turbu- 
lence with an integer order pole-zero TF fit and store the 
approximation in memory. 

In another embodiment of the present invention, a com- 
15 puter-implemented method includes simulating, by a proces- 
sor, atmospheric disturbances in both a time and frequency 
domain using an integer order pole-zero TF fit of a fractional 
order TF. The computer-implemented method also includes 
displaying, on a display device, a graphical representation of 
20 the integer order pole-zero TF fit. 

In yet another embodiment of the present invention, an 
apparatus includes physical memory including computer pro- 
gram instructions and a processor configured to execute the 
computer program instructions. The processor is configured 
25 to determine simulated atmospheric disturbances using an 
integer order pole-zero TF fit of a fractional order TF, and 
store the integer order pole-zero TF fit in the memory. 

BRIEF DESCRIPTION OF THE DRAWINGS 

30 

In order that the advantages of certain embodiments of the 
invention will be readily understood, a more particular 
description of the invention briefly described above will be 
rendered by reference to specific embodiments that are illus- 
35 trated in the appended drawings. While it should be under- 
stood that these drawings depict only typical embodiments of 
the invention and are not therefore to be considered to be 
limiting of its scope, the invention will be described and 
explained with additional specificity and detail through the 
40 use of the accompanying drawings, in which: 

FIG. 1 illustrates a graph of a longitudinal von Karman 
spectral and a final adjusted TF fit, according to an embodi- 
ment of the present invention. 

FIG. 2 illustrates a graph of a transverse von Karman 
45 spectral and a final adjusted TF fit, according to an embodi- 
ment of the present invention. 

FIG. 3 illustrates a graph of a temperature von Karman 
spectral and its TF fit, according to an embodiment of the 
present invention. 

50 FIG. 4 illustrates a graph of a von Karman acoustic wave 
velocity spectral due to temperature gust and its TF fits for 
different eddy dissipation rates (L=762 m), according to an 
embodiment of the present invention. 

FIG. 5 illustrates a graph of a von Karman acoustic wave 
55 velocity spectral due to temperature gust and its TF fits for 
different integral scale lengths (e=8.6e-5 m 2 /sec 3 ), according 
to an embodiment of the present invention. 

FIG. 6 illustrates a graph of a pressure von Karman spectral 
and its TF fit (€=8.6e-5 nr/sec 3 , L=762 m), according to an 
60 embodiment of the present invention. 

FIG. 7 illustrates a graph of a density von Karman spectral 
and its TF fit (e=8.6e-5 nr/sec 3 , L=762 m), according to an 
embodiment of the present invention. 

FIG. 8 illustrates a graph of a von Karman acoustic wave 
65 velocity spectral and its TF fits due to different values of eddy 
dissipation rates and integral length scales, according to an 
embodiment of the present invention. 
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FIG. 9 illustrates a graph of longitudinal acoustic wave 
velocity disturbance created using the TF fit of the longitudi- 
nal acoustic wave with unity input sinusoids, according to an 
embodiment of the present invention. 

FIG. 10 is a graph illustrating a thrust spectral density 5 
analysis, according to an embodiment of the present inven- 
tion. 

FIG. 11 is a flowchart illustrating a method of approximat- 
ing fractional order atmospheric turbulence with integer 
order pole-zero transfer function fits, according to an embodi- to 
ment of the present invention. 

FIG. 12 illustrates a computing system, according to an 
embodiment of the present invention. 

DETAILED DESCRIPTION OF THE is 

EMBODIMENTS 

Embodiments of the present invention simulate atmo- 
spheric disturbances in both the time and frequency domain 
using integer order TF fits of fractional order TFs. More 20 
specifically, such embodiments approximate fractional order 
atmospheric turbulence with integer order pole- zero transfer 
function fits for more accurate frequency and time domain 
simulations, von Karman type model forms of the Kolmog- 
orov spectral are constructed for the acoustic velocity gusts 25 
disturbances for both longitudinal and transverse distur- 
bances, and separately for ambient temperature and tempera- 
ture generated wind gusts, as well as pressure disturbances 
and density, if desirable. 

In some embodiments, a computer program that performs 30 
the modeling may run on an ordinary simulation development 
platform or any suitable compiler on a conventional com- 
puter. The simulation program may be configured to compute 
fractional order von Karman type TF approximations of the 
Kolmogorov spectra and compute the transfer function 35 
parameters associated with the atmospheric disturbance 
parameters, which may depend on the flight condition, 
vehicle speed, and selected atmospheric disturbance param- 
eters (e.g., worst case conditions). The simulation program 
may approximate the fractional order atmospheric turbulence 40 
with integer order pole-zero TFs, with the results displayed in 
the simulation environment in terms of the parameters of the 
various TFs and by plotting selected disturbances in the fre- 
quency domain. Also, the TFs may be simulated in the par- 
ticular simulation environment and appropriate time domain 45 
simulations may be developed, selected analyses may be 
conducted, and the results may be plotted. The flow logic of 
the computation process of the algorithm in such embodi- 
ments is shown in FIG. 11, which is discussed in more detail 
below. While different embodiments of the present invention 50 
can exist, with variations in the types of analyses and the 
results displayed, a feature of many embodiments is the meth- 
odology of scaling atmospheric disturbances into von Kar- 
man type forms and approximating their fractional order 
nature by computing products of integer order TFs to fit the 55 
fractional order. 

Because of the fractional order of these models, a circuit 
equivalent of atmospheric disturbances with associated 
capacitance and resistances, computed as a function of atmo- 
spheric disturbance parameters, is developed and used as the 60 
basis to derive the integer order pole-zero approximations. 
Utilizing such formulations, the transfer function poles and 
zeros of these approximations to the atmospheric distur- 
bances can be directly computed for different parameters 
describing atmospheric turbulence. The model is derived to 65 
duplicate the fractional order form of real atmospheric turbu- 
lence and is more accurate than conventional models. 


4 

The models may be computed and simulations may be 
performed over a wide range of altitudes and a wide degree of 
variations in atmospheric turbulence conditions. Some 
embodiments may be relatively easy to implement and may 
be highly representative of the actual fractional order nature 
of atmospheric turbulence. Embodiments may be applicable 
to propulsion system flow field type disturbances and vehicle 
gust loads. These new models could be used to design con- 
trols for aerospace vehicle propulsion or flight control sys- 
tems. 

Fractional Order Fit of Atmospheric Turbulence Model 

Atmospheric disturbance has a fractional order in nature, 
meaning that rather than a whole number, the order is a 
fraction such as 5/3. Fractional order differential equations 
are difficult to explicitly ornumerically solve. This is because 
fractional order derivatives, unlike integer order derivatives, 
do not obey the locality law (i.e., the limit theorem). Rather, 
fractional order derivatives obey the law of non-locality. This 
property makes the state transition matrix a convolution inte- 
gral, with each integration step necessitating the computation 
of the convolution of the state all the way back to time zero 
where the state had an initial value of zero. 

Accordingly, embodiments of the present invention derive 
integer order pole-zero product TF approximations to simu- 
late fractional order atmospheric disturbances. The method- 
ology of embodiments of the present invention is derived such 
that integer order TF approximations of atmospheric distur- 
bance are symmetrically centered about the fractional order 
TF by interleaving integer order poles and zeros. Correction 
factors to the transfer function parameters have been devel- 
oped for the different type of disturbances so that a good fit is 
achieved for a density of one pole-zero pair per frequency 
decade. If for some reason a higher accuracy is desired, which 
is likely unnecessary for many implementations, the same 
process can be followed for increased pole-zero densities. In 
an embodiment (as shown in the figures), plots of the von 
Karman spectra and the respective TF approximations are 
directly generated to compare these accuracies. 

The idea behind the staircase (pole-zero) approximation of 
the fractional order is to increase the accuracy of the staircase 
approximation by increasing the density of pole-zero pairs 
per frequency decade so that with increased pole-zero den- 
sity, the steps become shorter, eventually collapsing to the 
straight line of the fractional order TF as the number of poles 
and zeros of the approximation increase to infinity. For such 
an approach to work, the frequency of the poles and zeros 
need to be related to the atmospheric disturbance parameters 
and also be derived in a way such that the staircase TF 
approximation is symmetrically centered about the fractional 
order TF. 

The Kolmogorov spectral represents the spectral density of 
a structured random field of atmospheric turbulence as 

SXi)=a,e 2/3 r 5/3 (1) 

where S f (k) signifies the type of disturbance, a, is a con- 
stant for each type of disturbance, € represents the eddy 
dissipation rate (m 2 /sec 3 ), and k represents the wave number 
(cycles/m). More specifically, cq=0.15 for the longitudinal 
wind velocity gust in the direction of flight (dimensionless), 
a. =0.2 for the vertical or horizontal wind velocity gust (di- 
mensionless), a r =0.39 for temperature disturbance 
(°K 2 s 2 m -2 ), and a r =5xl0 -4 (P/T) 2 for pressure disturbance 
(Pa 2 s 2 m -2 ), where P and T are in total atmospheric quantities. 

For an in-flight vehicle encountering atmospheric turbu- 
lence with a wavenumber k (in cycles/m), the disturbance 
frequency (up to) experienced by the vehicle may be deter- 
mined by 
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where M is the vehicle Mach number, a is the local speed 
of sound (m/sec), which is altitude-dependent, and X (k=l/X) 
is the wavelength of the atmospheric turbulence (m/cycle). 
Later, Tank scaled the von Kannan spectrum to fit the Kol- 
mogorov model approaching the limit (i.e., for large values of 10 
k), which compares well with other data. For longitudinal 
distance, this spectrum is 


S LVK (k) = 2.7 e 2/3 L 5/3 


2 

[1 + (1.339(2tt)1/c) 2 ] 5 ^ 6 


(3) 15 


where k m „ is the adjustment factor, which equals 1 , but is 
represented in this form in case any adjustments need to be 
made to the atmospheric disturbance natural frequency. 

Utilizing these circuit analog parameters, a time domain 
disturbance can be derived based on an integer order TF 
approximation for atmospheric turbulence given the atmo- 
spheric disturbance parameters e, L, q, t, and r for the units 
conversion factor (with q=xr; x=5/3 and r =1 /3 for acoustic 
disturbances and l h for temperature and pressure) as 


P[ (*/0Jzi+ 1) 

= K uf „ -! w, 

m P 

n (s/io pi + 1 ) 

i 


( 8 ) 


For transverse disturbances, this spectrum is 


S,.vK(k) = 2.7 e 2/3 i 5/3 


1 + 5(1.339(2*)LU 2 
[1 + (1.339(2?r)Zjt) 2 ] 11/6 


20 


(4) 


In a nutshell, embodiments of the present invention were 
developed to fit the fractional order equations (3) and (4) with 
a product of first order poles and zeros symmetrically located 
on top of the fractional order for the acoustic waves. For 
temperature and pressure, equation (3) for the acoustic wave 3< 1 
is scaled by the differences in magnitude between the Kol- 
mogorov acoustic wave and the Kolmogorov temperature and 
pressure, respectively (see equation (1)). Then, the same 
approach is followed to fit the temperature and pressure spec- 
tral with products of first orders TFs. Finally, the fit TFs are 
adj usted with proportionality factors to improve the fit for the 
final TFs. Equation (2) is utilized to calculate the frequency 
based on the speed of the vehicle and size of the atmospheric 
turbulence, up to which unit amplitude sinusoids are con- 40 
structed as inputs to the TFs to provide for time domain 
simulations. 

The following description pertains to algorithms devel- 
oped with respect to embodiments of the present invention. 
The values of the equivalent capacitance and resistance of 45 
atmospheric turbulences may be determined by 


where W, is a series of sinusoids with unit amplitude fre- 
quency components distributed over the frequency range of 
interest, m_ is the number of zeros in the TF approximation, 
m p is the number of poles in the TF approximation, s is the 
Laplace operator, is the frequency of TF zero i of the 
integer order approximation (rad/sec), and (a pi is the fre- 
quency of TF pole of integer order approximation (rad/sec). 
The frequencies of the poles can be computed as 


w (9) 

Rwpi^Hpi j | (^Hpi I Mpi—j "it) 

7=i 

<u pi - 

10 * 2 '-‘>*«] ] (u Hpi / CO zH + 1 ) - 1 
7=1 

i = 2, 3 ... , m p 


where r| is the ratio of each decade interval where a pole or 
a zero will be used to estimate the fractional order TF and q is 
the fractional order of the equivalent electrical circuit, and 
where the first pole is computed as follows: 


« P i = ^(10* - 1) 


1 -? 


( 10 ) 


and the frequencies of the zeros can be computed as 


1 

( a,£ 2/i ) l,x (2nMa ) 


R,=1.3390(2x)(a/i 2,i ) Vx L (6) 

where C, is the equivalent electric circuit capacitance (far- 
ads) fort-type disturbance, R r is the equivalent electric circuit 
resistance (ohms) for t-type disturbance, a r is a constant for 
t-type disturbance, a is the speed of sound, 6 is the eddy 
dissipation rate (m 2 /sec 3 ), x is the fractional order of the 
atmospheric disturbance spectral, M is the Mach number, and 
L is the integral length scale (m). The corresponding natural 
frequency <n„ may be computed as 


7=1 

Wzi = : 

lO-^n (OlHa/hlp, + 1) - 1 

7=1 

i — 1,2 ... m z 

For a number n of frequency decades to be estimated, the 
number of poles and zeros in the TF approximation can be 
determined by 

rn p =(n- l)/r| (12) 


(5) 


50 


^ . n , , , , . 


(ID 


0) n 


Kum 


R,C, 


(7) 


m=m p -\ (13) 

For a desirable zero-pole density pair p pz per decade to be 
used to approximate the fractional order disturbance 
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n= 5 — 

2 Pp Z 

and the terms m Hp , and (o Hzi in equations (9) and (11) can be 
computed as 


=<a n (10' ng(2i ~ l ^-l) 1/q i=2 t 3 t . . . ,m p 

(15) 

=(O„(10 2T '« i -l) 1,9 i=2,3, . . . ,m z 

(16) 10 


The utility of u) Hpi and to Hzi is to maintain symmetry so that 
the staircase pole-zero approximation is symmetrically 
located on top of the fractional order TF asymptote. The 
proportional gains K OJ and K roz! in equations (9) and (11) are 
reserved for final adjustments that may be needed to these 15 
frequencies for more closely approximating the fractional 
order TFs representing atmospheric disturbance. 

Longitudinal and Transverse Acoustic Wave Turbulence 
The longitudinal and transverse acoustic wave atmospheric 
disturbances are in the form of pure wind gusts in the axial 211 
direction of vehicle motion and in the vertical direction of 
motion, respectively. Utilizing equations (5) through (16) 
with q=xr= s /9 (r=Vi), then in equation (8), based on 
inspection of the Tank model from equations (3) and (4), is 

^,=(5.4e 2/3 I 5/3 ) 1/3 (17) 

for the longitudinal disturbance, and 

iC ViS ,=(2.7€ 2/3 I 5/3 ) i/3 (18) 

for the transverse disturbance. For n=3 (i.e., the TF fit over 30 
the span of 3 decades in frequency), the proportionality factor 
adjustments to improve these fits have been found to be 

Ki,„=[K^„;K^ pl ;K^ t ]=[2A;\ 1 1/2.4 1/1.5;1 1 1] (19) 

for the longitudinal acoustic wave, and 35 

Jfv, u ,=/X™,Wio f >iWiazi]=[4.27;l 1 1/2.4 1/1.5;1 1 1] (20) 

for the transverse. These TF fits may be compared against 
the Tank von Karman spectral equations (3) and (4) with units 
conversion exponent r, shown in graphs 100 and 200 of FIGS. 40 
1 and 2 for e=8.6e-5 (m 2 /sec 3 ) and L=762 m. As shown in 
these figures, the fits provide a good approximation of these 
disturbances. More accuracy can be achieved around knees 
110 and 210 of these spectrals by increasing the density of 
pole-zero pairs per frequency decade in this region. While 45 
there may be a certain level of inaccuracy at lower frequen- 
cies, this is acceptable since typical control system designs 
have no problem attenuating disturbances at low frequency 
ranges (i.e., for the frequencies shown in the figures around 
the knees and below). 50 

Temperature Turbulence 

Atmospheric temperature turbulence causes both tempera- 
ture and acoustic velocity disturbances due to changes in the 
speed of sound. For the most part, temperature disturbances 
result in vertical displacement of air parcels (i.e., “gravity 55 
waves”). Thus, acoustic velocity disturbances caused by tem- 
perature will generate vertical wind gusts that add with any 
transverse acoustic velocity gusts. For a propulsion system, 
the wing forebody will turn a vertical gust into a longitudinal 
gust by multiplying the vertical gust by the conversion factor 60 


(M„ 




65 

This factor amounts to 0.63 for M CCi =2.35, and can be ignored 
for worst-case purposes. 


8 

Temperature turbulence spectral density also follows the 
5/3 law, but its units are in terms of Kelvin squared. Thus, in 
order to convert the units to Kelvin, the exponent r becomes 
V2, which makes the fractional order q=%. Based on this, the 
TF fit performed for temperature is done in a similar fashion 
as with the acoustic waves discussed in the previous section 
(i.e., utilizing equations (5) through (16)), with additional 
relations. 

K TfiK ,e mp r^ 14.0e 2,3 £ 5 ' 3 (21) 

1 1/1.1 1/1. 2;1 1 1] (22) 

These relations have been derived by utilizing the Kolmog- 
orov spectral of equation (1) to compute the difference in 
magnitude between the Kolmogorov acoustic wave and that 
of the temperature, then scaling the von Karman type tank 
model of equation (3) by this difference in magnitude in order 
to generate the von Karman type model for temperature. 

Graph 300 of FIG. 3 shows this TF fit with the Tank spectral 
of equation (3) (i.e., equation (3) substituting the scaled equa- 
tion (21) for the numerator and by also applying the units 
conversion factor in the denominator for e=8.6e-5 (m 2 /sec 3 ) 
and L=762 m). The proportionality factors K, for all of the fits 
can also be considered part of the formulations. Their accu- 
racy at higher frequencies starts deviating somewhat with 
higher values of L (integral length scales of a few thousand 
meters and above), values that are considered to be outside of 
the typical range. L=762 m is considered to be standard in the 
aerospace industry. 

Temperature turbulence also causes acoustic wave gusts 
due to change of the local speed of sound. The speed of sound 
is related with temperature as 

a^yRT, (23) 

where T z is the static temperature in degrees Kelvin and y is 
the ratio of specific heats (-1 .41). By perturbing the relation 
of the speed of sound and temperature in equation (23) and 
substituting the resulting expression for Aa into the Mach 
number equation of M=v/a, the relation between the change 
in temperature and acoustic velocity can be obtained as 


Av = 


MyR 

2 aT 


AT 


(24) 


This relation can be applied to equation (21) to compute 
^Tfit(acoustic) f° r the acoustic wave velocity disturbance due to 
an atmospheric temperature fluctuation as 


MyR 1 

KT.Macoumc) = 


(25) 


Equation (25) is applicable to transverse or vertical wind 
gust disturbances generated by temperature turbulence, and 
accurate in the worst case sense for longitudinal wind gusts, 
especially at higher supersonic Mach numbers (i.e., for Mach 
numbers higher than 2.35, the wing forebody in a propulsion 
system, mentioned earlier as the conversion factor, turns a 
transverse to a longitudinal disturbance using a factor calcu- 
lated to be 0.63 at Mach 2.35). Generally, longitudinal wind 
gusts due to temperature variations can be expressed as fol- 
lows by using this conversion factor 
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- 1 )yR i — — 7 — 

KT.fttUongitudinal) = . V 14.0 fi 2 / 3 /. 5 / 3 

2a oy jMl-l 


(26) 5 


This TF fit for the acoustic velocity disturbance at super- 
sonic cruise altitude, at Mach 2.35, for L=762 m and two 
different values of e (5e-4 and 8.6e-5) is shown in graph 400 
of FIG. 4. TheTF fit fore=8.6e-5 and for two different values 
of integral scale lengths is shown in graph 500 of FIG. 5. As 
can be seen from FIGS. 4 and 5 compared to FIG. 1, the 
acoustic velocity gusts produced by fluctuations in tempera- 1 . 
ture are much higher in amplitude than those produced by 
pure acoustic velocity gusts. 

For the combined acoustic wave, the acoustic wave veloc- 
ity due to temperature gust dominates in the low frequency 
range and at higher frequencies both acoustic waves affect the 
total. In the worst case, it can be assumed that both effects of 
the acoustic wave velocity can combine. It is also assumed 
that both the longitudinal and transverse acoustic wave 
velocities will not combine to produce worst case affects. For 
engines under the wing, the transverse wave can be assumed 
to be converted to a longitudinal wave via the wing forebody. 

Pressure Turbulence 

To develop the von Karman type model for pressure, the 
same scaling process as described in the temperature section 
may be followed for pressure. The units conversion exponent , f) 
for pressure is the same as that for temperature (i.e., r=Vi ), 
which makes the fractional order the same (i.e., q= 5 /6). As 
such, K P jj t and K P co may be computed as follows 

Kpjjpl 1 1 ,6e 2/3 Z 5/3 (27) 

35 

Kp^/K^K^K^M 1.5;1 1 1/1.1 1/1. 2;1 1 1] (28) 

Graph 600 of FIG. 6 shows this TF fit fore=8.6e-5 m 2 /sec 3 
and L=762 m compared to the scaled von Kannan type spec- 
tral (i.e., equation (27) substituted for the numerator and 
proportional factor of equation (3) and the denominator of 
equation (3), raised to the power r). 

Density Turbulence 

Density disturbances are contained within temperature and 
pressure fluctuations. As such, its inclusion in a simulation, 
along with pressure and temperature, will produce unneces- 
sary additional disturbances. However, if need be, the relation 
of density disturbance with temperature and pressure can be 
derived by perturbing the state equation as 


PpAP Pq 
R(T 0 + AT) RTo 


(29) 


A plot of density disturbance utilizing equation (29) is 55 
shown in graph 700 of FIG. 7. 

Simplified Atmospheric Turbulence Model for Controls 
Design 

The fractional order TF fit may be further simplified by 
identifying certain characteristics of this model that are 60 
important for controls design in order to handle atmospheric 
turbulence. Based on observations utilizing FIG. 4, differ- 
ences in eddy dissipation rates e produce spectral densities 
with the same frequencies, offset by a fixed magnitude. This 
offset is accounted by differences in TF gain. 65 

Calculating these frequencies using the formulations 
developed so far also supports this observation. Differences 
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in the fractional order of the disturbances, such as longitudi- 
nal acoustic wave velocity versus temperature will produce 
different spectral frequencies. Unlike differences in eddy dis- 
sipation rates that produce the same frequencies with differ- 
ent TF magnitudes, inspection of FIG. 5 shows that spectral 
with different integral length scales L may be represented by 
different frequencies. However, different integral length 
scales may only affect the low frequency spectral, where a 
typical feedback controls design should have no difficulty 
attenuating disturbances at such a low frequency range. 

Graph 800 of FIG. 8 shows spectral density TF fits, first, 
with a fixed integral length scale and two different values of E. 
This shows that the differences in e produce a fixed offset in 
magnitude, but with the same pole-zero frequencies. FIG. 8 
also shows a comparison between two spectral densities with 
the same value of e and two different values of L. This com- 
parison, with different values of L, is represented with differ- 
ent pole-zero frequencies. However, the high frequency 
asymptotes match, which means that for the purpose of con- 
trols design, either of these TF fits can be utilized despite their 
differences in pole-zero frequencies. 

Therefore, the frequency domain atmospheric turbulence 
model can be simplified by utilizing the model developed in 
the previous section to calculate the fixed frequencies due to 
a fixed integral length scale and variable eddy dissipation 
rates. Based on that, the simplified model developed for atmo- 
spheric turbulence, which ignores differences in the pole-zero 
frequencies due to different integral length scales can be 
represented as 


G LA (s) = 70sU 9 - 
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Wliere G LA , G VA , G^, G TLA , G TVA , and G P are the simplified 
atmospheric disturbance TFs for longitudinal and transverse 
acoustic wave velocity, temperature, longitudinal acoustic 
velocity due to temperature, and vertical acoustic velocity 
due to temperature, and pressure, respectively. The factor 

(M„ - 1) / J Ml - 1 

in equation (33) can be set for worst case turbulence condi- 
tions, especially for relatively high Mach numbers, in which 
case it makes equation (33) and equation (34) the same. It is 
noted that M„ pertains to the freestream Mach number, which 
is the same as the vehicle Mach number. 

The TF fit described in the previous section was carried up 
to 200 FIz, which should make it accurate even beyond this 
frequency range. However, the model simplification repre- 
sented by equations (30) to (35) introduces additional errors 
at high frequencies and its validity should be limited to about 
200 Hz. The frequencies of this model are based on an integral 
length scale of L=762 m. This means that this simplified 
model may not be unique, and other equivalent models can be 
calculated using different values of L with the equations 
developed in the previous sections. Thus, for a controls design 
which can adequately attenuate 1 Hz disturbances, this model 
is accurate for values of L as low as 100 m, which is well 
below the expected range of L. A value of L=762 m is con- 
sidered standard in the aerospace industry, as mentioned pre- 
viously. 

Time Domain Disturbances 

The TF fits in the previous section can be thought of as 
filters, which take input disturbances at desirable frequencies 
and convert them to representative free stream atmospheric 
disturbances. The TF fits contain the representative magni- 
tudes of these disturbances. Thus, the discrete time domain 
frequencies W r may only need to have unity root mean square 
(RMS) values. For instance, a time domain input can be 
constructed as a sum of unit amplitude sinusoids, starting at 
the low frequency asymptote of the spectral, and continuing 
up with sinusoids that are dispersed along the high frequency 
asymptote up to frequencies that cover the controller band- 
width. Also, the control design can be checked with single 
frequency sinusoids, sequentially simulated. The total energy 
or the maximum amplitude of the disturbance will need to be 
reasonable. For instance, it may be considered that total wind 
gust will not exceed an amplitude of 1 80 mph, or some other 
reasonable number that depends on altitude. 

For most of the spectral density simulations covered so far, 
the worst case eddy dissipation rate for low to moderate 
turbulence was used for altitudes of about 36,000 ft (~1 1 km) 
and above. However, the highest atmospheric turbulence ever 
recorded had an eddy dissipation rate of 1 .7e-3 (m 2 /sec 3 ). 
Theoretically, it is possible that even more severe turbulence 
can be encountered. 

Conventional approaches cover atmospheric turbulence 
probabilistically, as an exceedance for controls design pur- 
poses. In terms of probabilistic atmospheric turbulence mod- 
els, there are no comparisons drawn here between those mod- 
els and the non-probabilistic approach described herein, even 
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though exceedance can also be incorporated with this 
approach. The modeling approach developed herein is geared 
towards implementing a controls design by considering worst 
case atmospheric turbulence. This also allows implementa- 
5 tion of control designs to cover worst case turbulence encoun- 
ters without penalizing performance under normal atmo- 
spheric conditions by incorporating control scheduling in the 
design. 

Graph 900 of FIG. 9 shows a time domain longitudinal 
to acoustic wave disturbance utilizing the acoustic wave TF fit 
discussed in the previous section. The input to this TF is a 
combination of unity amplitude sinusoids distributed along 
the frequency range shown in FIG. 1. For a single unity 
amplitude sinusoid at low frequency, the amplitude of the 
15 acoustic disturbance would be approximately 9 m/sec accord- 
ing to FIG. 1. For this summation of unity amplitude input 
sinusoids, their amplitudes combine to give approximately 15 
msec maximum disturbance. The same approach can be used 
to construct other atmospheric disturbances for time domain 
20 simulations. 

FIG. 1 0 is a graph 1000 illustrating a thrust spectral density 
analysis, according to an embodiment of the present inven- 
tion. The thrust spectral density analysis is for a certain pro- 
pulsion system subjected to atmospheric disturbances fortwo 
25 different eddy dissipation rates (severity of disturbances) 
with unity frequency sinusoids selected as inputs to the atmo- 
spheric disturbance TFs to limit the combined wind gust 
velocities to certain values. Graph 1000 also shows the spec- 
tra density limits with frequency, with the dashed line indi- 
30 eating the worst anticipated bounds of tins spectral. Even 
though not shown in FIG. 10, a different choice of input 
sinusoids of lower frequencies that also meets the selected 
wind gust limits will produce spectral densities with ampli- 
tudes with peaks that are bounded by these straight line 
35 bounds at these lower frequencies. Such an analysis, for 
instance, can be used to assess the vehicle ride quality or to 
assess the propulsive forces coupled to the body modes of the 
vehicle. 

FIG. 11 is a flowchart 1100 illustrating a method of 
40 approximating fractional order atmospheric turbulence with 
integer order pole-zero transfer function fits, according to an 
embodiment of the present invention. In some embodiments, 
the method of FIG. 11 may be implemented, for example, by 
computing system 1200 of FIG. 12. The method begins with 
45 selecting the flight conditions and disturbance parameters at 
1105. Next, the pole- zero density and number of frequency 
decades are selected at 1110. The corresponding natural fre- 
quency is computed at 1 1 15 and the pole-zero frequencies are 
computed at 1120. 

50 If simplified TFs are to be used at 1125, simplified TFs for 
controls design are computed at 1130. Otherwise, adjusting 
factors are applied and TFs are computed for each disturbance 
type at 1135. The TFs are then configured for the overall 
disturbance simulation at 1140, and unity amplitude sinuso- 
55 ids are constructed and distributed up to a frequency corre- 
sponding to the vehicle speed at 1145. The simulation is then 
run at 1150. 

If the selected maximum wind velocity is not satisfied at 
1155, the process returns to adjusting the unity amplitude 
60 sinusoids (i.e., appropriately changing the frequencies of the 
sinusoids and/or the number of sinusoids) at 1145. Otherwise, 
the process checks whether the user desires to generate 
another model at 1160. If so, the process returns to the begin- 
ning at 1105. Otherwise, the model is run and the results are 
65 plotted at 1165. Optionally, cases may be run and spectral 
density analysis may be constructed at 1170. Otherwise, the 
process ends or the simulation is connected to the vehicle 
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system model to run simulations. In some embodiments, the 
algorithm process is setup to run automatically in a computer 
program, and some of the steps like setting up the process, 
choosing the unit amplitude sinusoids, generating a different 
model, and running vehicle system analysis may require user 5 
involvement. 

FIG. 12 illustrates a computing system, according to an 
embodiment of the present invention. System 1200 includes a 
bus 1205 or other communication mechanism for communi- 
cating information, and a processor 1210 coupled to bus 1205 to 
for processing information. Processor 1210 may be any type 
of general or specific purpose processor, including a central 
processing unit (CPU) or application specific integrated cir- 
cuit (ASIC). System 1200 further includes a memory 1215 for 
storing information and instructions to be executed by pro- 15 
cessor 1210. Memory 1215 can be comprised of any combi- 
nation of random access memory (RAM), read only memory 
(ROM), flash memory, cache, static storage such as a mag- 
netic or optical disk, or any other types of non-transitory 
computer-readable media or combinations thereof. Addition- 20 
ally, system 1200 includes a communication device 1220, 
such as a wireless network interface card, to provide access to 
a network. 

Non-transitoty computer-readable media may be any 
available media that can be accessed by processor 1210 and 25 
may include both volatile and non-volatile media, removable 
and non-removable media, and communication media. Com- 
munication media may include computer-readable instruc- 
tions, data structures, program modules or other data in a 
modulated data signal such as a carrier wave or other transport 30 
mechanism and includes any information delivery media. 

Processor 1210 is further coupled via bus 1205 to a display 
1225, such as a Liquid Crystal Display (“LCD”), for display- 
ing information to a user. A keyboard 1230 and a cursor 
control device 1235, such as a computer mouse, are further 35 
coupled to bus 1205 to enable a user to interface with system 
1200. However, in certain embodiments such as those for 
mobile computing implementations, a physical keyboard and 
mouse may not be present, and the user may interact with the 
device solely through display 1225 and/or a touchpad (not 40 
shown). Any type and combination of input devices may be 
used as a matter of design choice. 

In one embodiment, memory 1215 stores software mod- 
ules that provide functionality when executed by processor 
1210. The modules include an operating system 1240 for 45 
system 1200. The modules further include an atmospheric 
disturbance simulation module 1245 that is configured to 
simulate atmospheric disturbances using a TF fit. System 
1200 may include one ormore additional functional modules 
1250 that include additional functionality. 50 

One skilled in the art will appreciate that a “system” could 
be embodied as a personal computer, a server, a console, a 
personal digital assistant (“PDA”), a cell phone, a tablet com- 
puting device, or any other suitable computing device, or 
combination of devices. Presenting the above-described 55 
functions as being performed by a “system” is not intended to 
limit the scope of the present invention in any way, but is 
intended to provide one example of many embodiments of the 
present invention. Indeed, methods, systems and apparatuses 
disclosed herein may be implemented in localized and dis- 60 
tributed fonns consistent with computing technology, includ- 
ing cloud computing systems. 

It should be noted that some of the system features 
described in this specification have been presented as mod- 
ules, in order to more particularly emphasize their implemen- 65 
tation independence. For example, a module may be imple- 
mented as a hardware circuit comprising custom very large 
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scale integration (“VLSI”) circuits or gate arrays, off-the- 
shelf semiconductors such as logic chips, transistors, or other 
discrete components. A module may also be implemented in 
programmable hardware devices such as field programmable 
gate arrays, programmable array logic, programmable logic 
devices, graphics processing units, or the like. 

A module may also be at least partially implemented in 
software for execution by various types of processors. An 
identified unit of executable code may, for instance, comprise 
one or more physical or logical blocks of computer instruc- 
tions that may, for instance, be organized as an object, proce- 
dure, or function. Nevertheless, the executables of an identi- 
fied module need not be physically located together, but may 
comprise disparate instructions stored in different locations 
which, when joined logically together, comprise the module 
and achieve the stated purpose for the module. Further, mod- 
ules may be stored on a computer-readable medium, which 
may be, for instance, a hard disk drive, flash device, random 
access memory (“RAM”), tape, or any other such medium 
used to store data. 

Indeed, a module of executable code could be a single 
instruction, or many instructions, and may even be distributed 
over several different code segments, among different pro- 
grams, and across several memory devices. Similarly, opera- 
tional data may be identified and illustrated herein within 
modules, and may be embodied in any suitable form and 
organized within any suitable type of data structure. The 
operational data may be collected as a single data set, or may 
be distributed over different locations including over different 
storage devices, and may exist, at least partially, merely as 
electronic signals on a system or network. 

The method steps performed in FIG. 1 1 may be performed 
by a computer program, encoding instructions for the nonlin- 
ear adaptive processor to perform at least the method 
described in FIG. 11, in accordance with an embodiment of 
the present invention. The computer program may be embod- 
ied on a non-transitory computer-readable medium. The com- 
puter-readable medium may be, but is not limited to, a hard 
disk drive, a flash device, a random access memory, a tape, or 
any other such medium used to store data. The computer 
program may include encoded instructions for controlling the 
nonlinear adaptive processor to implement the method 
described in FIG. 11, which may also be stored on the com- 
puter-readable medium. 

The computer program can be implemented in hardware, 
software, ora hybrid implementation. The computer program 
can be composed of modules that are in operative communi- 
cation with one another, and which are designed to pass 
information or instructions to display. The computer program 
can be configured to operate on a general purpose computer, 
or an application specific integrated circuit (“ASIC”). 

Some embodiments of the present invention use deriva- 
tions to approximate fractional order atmospheric turbulence 
with integer order pole-zero transfer function fits for more 
accurate time domain simulations. An existing von Karman 
type model form of the Kolmogorov spectral is utilized ini- 
tially for acoustic velocity gust disturbances. Hiis model may 
be scaled in order to develop the von Karman forms for other 
types of atmospheric disturbance, such as temperature, pres- 
sure, and density. A circuit equivalent may be used as the basis 
to derive the integer order pole- zero approximations. The TF 
poles and zeros of the approximations to the atmospheric 
disturbances can be directly computed for different param- 
eters describing atmospheric disturbance. 

It will be readily understood that the components of various 
embodiments of the present invention, as generally described 
and illustrated in the figures herein, may be arranged and 
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designed in a wide variety of different configurations. Thus, 
the detailed description of the embodiments of the systems, 
apparatuses, methods, and computer programs of the present 
invention, as represented in the attached figures, is not 
intended to limit the scope of the invention as claimed, but is 5 
merely representative of selected embodiments of the inven- 
tion. 

The features, structures, or characteristics of the invention 
described throughout this specification may be combined in 
any suitable manner in one or more embodiments. For to 
example, reference throughout this specification to “certain 
embodiments,” “some embodiments,” or similar language 
means that a particular feature, structure, or characteristic 
described in connection with the embodiment is included in at 
least one embodiment of the present invention. Thus, appear- 15 
ances of the phrases “in certain embodiments,” “in some 
embodiment,” “in other embodiments,” or similar language 
throughout this specification do not necessarily all refer to the 
same group of embodiments and the described features, struc- 
tures, or characteristics may be combined in any suitable 20 
manner in one or more embodiments. 

It should be noted that reference throughout this specifica- 
tion to features, advantages, or similar language does not 
imply that all of the features and advantages that may be 
realized with the present invention should be or are in any 25 
single embodiment of the invention. Rather, language refer- 
ring to the features and advantages is understood to mean that 
a specific feature, advantage, or characteristic described in 
connection with an embodiment is included in at least one 
embodiment of the present invention. Thus, discussion of the 30 
features and advantages, and similar language, throughout 
this specification may, but do not necessarily, refer to the 
same embodiment. 

Furthermore, the described features, advantages, and char- 
acteristics of the invention may be combined in any suitable 35 
manner in one or more embodiments. One skilled in the 
relevant art will recognize that the invention can be practiced 
without one or more of the specific features or advantages of 
a particular embodiment. In other instances, additional fea- 
tures and advantages may be recognized in certain embodi- 40 
ments that may not be present in all embodiments of the 
invention. 

One having ordinary skill in the art will readily understand 
that the invention as discussed above may be practiced with 
steps in a different order, and/or with hardware elements in 45 
configurations which are different than those which are dis- 
closed. Therefore, although the invention has been described 
based upon these preferred embodiments, it would be appar- 
ent to those of skill in the art that certain modifications, 
variations, and alternative constructions would be apparent, 50 
while remaining within the spirit and scope of the invention. 

In order to determine the metes and bounds of the invention, 
therefore, reference should be made to the appended claims. 

NOMENCLATURE 55 

The following Nomenclature is used in claim 6 and is defined 
as: 

q fractional order of equivalent electrical circuit 

r| the ratio of each decade interval where a pole or a zero will 60 

be used to estimate the fractional order TF 
rm number of zeros in the TF approximation 
m p number of poles in the TF approximation 
( » . frequency of TF pole i of integer order approx, (rad/sec) 
u) r! frequency of TF zero i of integer order approx, (rad/sec) 65 
K (0/)I proportional gain for the calculation of the ith pole 

frequency 


(x) Hpl symmetry frequency of the ith zero (rad/sec) to guaran- 
tee symmetric placement of approximation 

proportional gain for the calculation of the ith zero fre- 
quency 

a ) Hzi symmetry frequency of the ith zero (rad/sec) to guaran- 
tee symmetric placement of approximation. 

The invention claimed is: 

1. A non-transitory computer readable storage storing 
computer executable instructions which when executed by a 
computer perform a process comprising: 

scaling a Kolmogorov spectral to convert the Kolmogorov 
spectral into a finite energy von Karman spectral; 
deriving a fractional order pole-zero transfer function (TF) 
from the von Karman spectral; 
approximating fractional order atmospheric turbulence 
with an integer order pole-zero TF fit; and 
storing the approximation in memory. 

2. The process of claim 1, further comprising: 
displaying, via a display device, a graph of the integer order 

pole- zero TF fit of the von Karman spectral. 

3. The process of claim 1, further comprising: 
receiving and incorporating a user-selected flight condition 

and atmospheric turbulence parameters. 

4. The process of claim 1, wherein the integer order pole- 
zero TF fit is symmetrically centered about the fractional 
order pole-zero TF by interleaving integer order poles and 
zeros. 

5. The process of claim 1, wherein a frequency of integer 
order poles and zeros is related to atmospheric disturbance 
parameters derived such that the integer order pole-zero TF fit 
is symmetrically centered about the fractional orderpole-zero 
TF. 

6 . The process of claim 1, wherein the process further 
causes the processor to determine a frequency of poles by an 
equation 


Kuipi £ OHpi | | {(JJHpi l ^ pi-j 4 1 ) 

M . . . 

oj pi = — i = 2,3 ... ,m p 

jOtfa-Wn (.<•>&/ o>d-j + D-l 
./=! 

and determine a frequency of zeros by an equation 


0J n-j + 1) 

M ■ o 

co p = : i = 2, 3 ... , m z . 

10- 2 ^n («W*V + 1)-1 

M 

wherein 

q=fractional order of equivalent electrical circuit, 
r|=the ratio of each decade interval where a pole or a zero 
will be used to estimate the fractional order TF, 
rrL=number of zeros in the TF approximation, 
m p number of poles in the TF approximation, 

^frequency of TF pole i of integer order approx, (rad / 
sec), 

<B ZJ =frequency of TF zero i of integer order approx, (rad / 
sec), 

K(opi=proportional gain for the calculation of the ith pole 
frequency, 

dissymmetry frequency of the ith zero (rad/sec) to guar- 
antee symmetric placement of approximation, 
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K a)e! =proportional gain for the calculation of the ith zero 
frequency, 

to^=symmetry frequency of the ith zero (rad/ sec) to guar- 
antee symmetric placement of approximation. 

7. The process of claim 1, wherein the process further 5 
causes the processor to determine a longitudinal fit by equa- 
tion 

^ r (5.4e 2,3 Z 5,3 ) 1/3 

determine a transverse disturbance by to 

'S: v ^=(2.7e 2,3 Z. 5,3 ) 1/3 , 

and scale associated TF frequencies for a better fit, wherein 
e=eddy dissipation rate (m 2 /sec 3 ), 

L=integral length scale. 15 

8. The process of claim 1, wherein the process further 
causes the processor to determine an integer order pole-zero 
TF fit for temperature by equation 

K^pr'Iu.O^L™ 20 

wherein 

e=eddy dissipation rate (m 2 /sec 3 ), 

L=integral length scale, 

scale associated integer order pole-zero TF parameters for 
a better fit, and 

detennine the integer order pole-zero TF parameters asso- 
ciated with wind gusts due to temperature and pressure 
disturbances. 

9. A computer-implemented method, comprising: 
simulating, by a processor, atmospheric disturbances in 

both a time and frequency domain using an integer order 
pole- zero transfer function (TF) fit of a fractional order 
pole- zero TF; and 

displaying, on a display device, a graphical representation 
of the integer order pole-zero TF fit. 

10. The computer-implemented method of claim 9, further 
comprising: 

deriving, by the processor, the fractional order pole-zero 
TF from a von Karman spectral. 

11. The computer-implemented method of claim 10, 
wherein the von Karman spectral is a finite energy conversion 
of a Kolmogorov spectral. 

12. The computer-implemented method of claim 10, fur- 
ther comprising: 
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selecting, by the processor, an offset decibel level such that 
a resultant horizontal asymptote of the integer order 
pole-zero TF fit approximates a corresponding horizon- 
tal asymptote of the von Karman spectral. 

13. The computer-implemented method of claim 9, 
wherein the integer order pole-zero TF fit is symmetrically 
centered about the fractional order pole- zero TF by interleav- 
ing integer order poles and zeros. 

14. The computer-implemented method of claim 9, 
wherein a frequency of integer order poles and zeros is related 
to atmospheric disturbance parameters derived such that the 
integer order pole-zero TF fit is symmetrically centered about 
the fractional order pole- zero TF. 

15. An apparatus, comprising: 

physical memory comprising computer program instruc- 
tions; and 

a processor configured to execute the computer program 
instructions, the processor configured to: 
determine simulated atmospheric disturbances using an 
integer order pole-zero transfer function (TF) fit of a 
fractional order pole-zero TF, and 
store the integer order pole-zero TF fit in the memory. 

16. The apparatus of claim 15, the processor further con- 
figured to: 

derive the fractional order pole-zero TF from a von Karman 
spectral. 

17. The apparatus of claim 16, wherein the von Karman 
spectral is a finite energy conversion of a Kolmogorov spec- 
tral. 

18. The apparatus of claim 16, the processor further con- 
figured to: 

select an offset decibel level such that a resultant horizontal 
asymptote of the integer order pole-zero TF fit approxi- 
mates a corresponding horizontal asymptote of the von 
Karman spectral. 

19. The apparatus of claim 15, wherein the integer order 
pole-zero TF fit is symmetrically centered about the frac- 
tional order pole-zero TF by interleaving integer order poles 
and zeros. 

20. The apparatus of claim 15, wherein a frequency of 
integer order poles and zeros is related to atmospheric distur- 
bance parameters derived such that the integer order pole- 
zero TF fit is symmetrically centered about the fractional 
order pole-zero TF. 





